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FOREWORD 


The purpose of this report is to present a description of the 
computational fluid dynamics capability of the General Interpolants 
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Company, Huntsville Research and Engineering Center. It was 
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1. INTRODUCTION 


In recent years, extensive literature has appeared on the numerical 
integration of some form of the Navier -Stokes equations of fluid dynamics. 
The introduction of larger and faster computers has enabled investigators 
to experiment numerically with various approximate solution procedures 
for these complex equations. This development has proceeded along two 
different, but related, avenues: (1) mathematical solution of partial differ- 
ential equations, and (2) engineering computation of physical phenomena. 
Although item one is not ignored, we are primarily concerned here with 
engineering analysis for obtaining useful results in practice. Acceptance 
of unsteady techniques in conjunction with larger and faster computers is 
changing the philosophy of computational mechanics since it is now possible 
to begin thinking of a general prediction tool for computing real world prob- 
lems . 


The system of equations considered here are various subsets of the 
Navier -Stokes equations plus an appropriate energy equation. Lagrangian 
coordinates or mixed Eulerian- Lagrangian frames are not considered. The 
equations are cast in unsteady form to pose the problem as initial valued in 
time. With appropriate spatial boundary values, the equations can be inte- 
grated forward in time in a parabolic/hyperbolic sense. The use of the un- 
steady formulation has a number of advantages: (1) the transient portion of a 
flow can be obtained as well as a steady state; (2) the asymptotic approach 
to a steady flow proceeds in much the same manner as nature; (3) subsonic, 
transonic and supersonic flows can be calculated with same methodology; 

(4) multi -dimensional flows can be handled more conveniently; (5) with a 
proper form of the differential equations, integration can automatically 



proceed through shock waves or other discontinuities without undue compli- 
cation; and (6) computation of viscous flows with time -fluctuating turbulence 
can be performed more naturally. 

The art of computational mechanics is the construction of an appropriate 
numerical equivalent of the conservation equations since there is no systematic 
procedure which will yield the optimum approximation for a general system. 
There are at least two distinct ways to derive such approximation formulas. 

The first, or "grid" approach operates on the differential form of the equations 
by expanding the terms in a Taylor series and truncating it in some manner. 
The second approach, termed "cell" or "element" methods operate on the 
integral form of the equations and assumes a functional form for the varia- 
tion of the dependent variables over a local finite volume. In either approach, 
the concepts of a "consistent" approximation, an "accurate" calculation pro- 
cedure and a "stable" integration formula are all vital to obtaining solutions to 
the equations , 

The classical works use finite difference techniques (Ref. 1) but recently 
the concept of the finite element, borrowed from structural analysis (Ref. 2), 
has become increasingly popular. Scientists with extensive backgrounds in 
finite difference solutions defended their approach against the finite element 
advocates. Likewise the finite element proponents claimed that the other side 
did not understand their new method. The resulting arguments have produced 
many papers and few solutions. An occasional attempt has even been made to 
relate the two approaches, but the schism still remained. It is the intention of 
this discussion to produce a consistent approach which clarifies the situation, 
and as it turns out, demonstrates that both approaches have merit and that each 
camp can learn a great deal from the other. 

This report introduces a new methodology for constructing numerical 
analogs of the partial differential equations of continuum mechanics. A gen- 
eral formulation is provided which permits classical finite element methods 
and many of the finite difference methods to be derived directly. The approach, 
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termed the General Interpolants Method (GIM), is new in the sense that it can 
combine the best features of finite element and finite difference methods. The 
technique allows complex geometries to be handled in the finite element manner 
and operates on the integral form of the conservation laws. Solutions can be 
generated implicitly with the finite element analogs or by explicit finite differ- 
ence analogs, which do not require a reduction of large systems of linear alge- 
braic equations (no matrix inverse). A quasi-variational procedure is used to 
introduce boundary conditions into the method and to provide a natural assembly 
sequence for combining the element equations into the full domain equations. 

The domain of interest is first discretized by appropriate subdivision into 
an assemblage of interconnected finite elements. A mesh generation is used in 
the GIM approach which incorporates general curvilinear coordinates, stretching 
transformations and bivariate blending to produce an automated mesh/element 
generation. Shape functions based on a set of generalized interpolants are then 
chosen to describe the behavior over each element. We then proceed, as in the 
Method of Weighted Residuals (Ref. 3) by multiplying the discretized equations by 
a set of weight functions and integrating over the volume of the element. A quasi- 
variational procedure is then used to construct the assembled system of equations 
from the element equations, and to introduce boundary conditions into the method. 
By choosing the weight functions equal to the shape functions, we reproduce via 
Galerkin the classical finite element nodal analogs. It is at this point that we 
introduce one of the important concepts of GIM: orthogonal weight/shape functions. 
By appropriately choosing the weight functions to be orthogonal to the shape func- 
tions, we can obtain explicit nodal analogs. Further, by a coice of arbitrary con- 
stants in the orthogonal weight functions, we can reproduce known finite difference 
nodal analogs, such as centered difference, upwind/downwind differences and the 
two-step MacCormack algorithm. As a result of this spatial discretization, we 
have reduced the partial differential equations to ordinary differential equations 
with "time" as the independent variable. Any forward marching algorithm such 
as Euler, Runge-Kutta or predictor-corrector can be used to advance the solu- 
tion profiles in time. 

In this report, we will use the unsteady equations of two-dimensional in- 
viscid fluid flow to illustrate the GIM approach. Extension of this development 
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to three space dimensions and to viscous flow equations has been accomplished 
but we use the simpler problem to present the basic principles of GIM more 
tractably. The following points summarize GIM and its advantages over pre- 
vious approaches. 

• GIM is a higher order procedure for constructing compu- 
tational analogs of the conservation laws of continuum 
mechanics . 

• GIM is a total calculational procedure in that arbitrary 
geometries, expandable equation sets, different nodal 
analogs, general curvilinear coordinate systems and 
multi -constraint boundary values are handled. 

• Non-orthogonal sets of weight/shape functions give rise 
to analogs which have generally been classified as finite 
element methods. 

• Orthogonal sets of weight/shape functions give rise to 
explicit finite difference nodal analogs. 

• A single analysis can be generated which selectively 
employs either of several finite element approaches 

or either of several finite difference techniques accord- 
ing to the wishes of the analyst. 

• Finite Element discussions in the literature have mentioned 
qua si -variational methods, but none have consistently applied 
these ideas to boundary conditions as we have in GIM. 

• The classical finite element analogs, derived via Galerkin, 
are unconditionally unstable for solution of the strong con- 
servation form of the Euler equations for shock capturing. 

Several authors have presented formulations entitled Finite 
Element with automatic shock capturing, but no solutions 
are reported. We have found by numerical experiment, 
that the nodal analogs they derive are unstable. The only 
solutions presented make use of the theory of "weak solu- 
tions," which destroys the strong conservation form of the 
equations . 

• The GIM approach allows the flexibility and generality of 
finite element techniques to be effectively married with 
proven successful finite difference techniques to produce 
a superior higher order methodology. 
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2. METHOD DESCRIPTION 


2.1 EQUATION SETS AND THE QUASI- VARIATIONAL TECHNIQUE 

The partial differential equations governing the two-dimensional flow 
of an inviscid, nonconducting ideal gas can be written in the following popular 
form: 


9U 3E 

at + 3x 


8F 

ay 


+ H 


0 


( 1 ) 


where x, y are spatial Cartesian coordinates, t is time, U is a column vector 
of dependent flow variables, and E, F, H are nonlinear functions of U. For 
example , 
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where p = density, u = velocity component in x, v = velocity component in y, 
E = total energy per unit mass and P = pressure. 

This equation set, in conjunction with an equation of state P = P(p,E) 
and appropriate boundary values, formally defines the type of problem used 
here to describe GIM. 

This form of the equations is termed Strong Conservation Law or Di- 
vergence form. This form has been advantageous for finite difference pro- 
ponents because of the superior characteristics in the neighborhood of dis- 
continuities and steep gradients. It is also a convenient form for finite 
element analysis since the assembly procedure for collecting influence 
coefficients needs to be done only once. 
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GIM is developed in terms of equations of the form (1). Cartesian co- 
ordinates are used to write the equation sets and transformations are made 
numerically to relate a specific geometric region to this reference frame. 

Using Eq. (1) also allows the method to be formulated for multi -equation sets 
since the length of the column vectors is arbitrary. Thus additional equations 
can be added such as a third dimension, multi- specie reacting flow conservation 
laws and equations for multi-phase flows. Additional terms can also be added 
to this form to model viscous flows by including stress components, body forces, 
etc., to the E,F and H vectors. 

It is natural to view Eq. (1) as a time -dependent system. With proper 
initial values and a realistic time, this form indeed represents the transient 
behavior of the flow. Additional ways of viewing Eq. (1) are also possible. With 
arbitrary initial conditions, the time t can be considered only as a relaxation 
coordinate such that Eq. (1) represents an asymptotic steady state form. 

We may also view "t M as a spatial marching coordinate such that Eq. (1) 
represents a steady state hyperbolic system. Although a spatial marching varia- 
tion of GIM is not precluded, we consider here only the unsteady or relaxation 
viewpoint of Eq. (1). 

In many problems in mathematics and physics, a variational principle 
exists in which a stationary value of the principal function yields the governing 
equations. Perhaps the most complex, yet familiar, example is that of the maxi- 
mization of entropy for closed systems leading to the equations governing chem- 
ical equilibrium. 

There are few examples of field solutions obtainable through variational 
principles; the most notable in fluid mechanics is restricted to potential flow. 
There is no known principle that generates the primitive equations of motion 
for fluid mechanics. 

In spite of the absence of a known principle we may, however, still take 
advantage of variational techniques. Let us assume a functional exists: 


6 



I 



whose extremum yields the primitive equation of motion. That is 



At the extremum, therefore, the coefficient of the variations must vanish, 
thereby satisfying the equations of motion. The above statements are described 
as a quasi-variational technique and where, in the absence of the variational 
principle, provide a rigorous academic point of departure for further develop- 
ments. This type of idea is used in the GIM development to form an assembly 
procedure and to describe multi-constraint boundary values. 
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2.2 GEOMETRIC TREATMENT 


The domain of interest is considered to be geometrically arbitrary in that 
any shape is represented as a bivariate blend of regular subdomains. Figure 
2-1 shows an example of a full geometric domain in which flow is to be computed. 
This domain is subdivided into four regions as shown. In general this division 
is made such that analytical functions describe the shape of each side but it can, 
of course, be made by point specification and piecewise linear sides. Attention 
can then be focused on each region with global coordinates (x, y) and local region 
coordinates (£’,1T). The regions are blended together at the junctions to provide 
the continuous full domain geometry. 

We then proceed by considering each region separately as depicted in 
Fig. 2-2. Gordon and Hall (Ref. 5) show that the general relationship between 
physical space and local curvilinear coordinates is given by the following 
transformation: 

R (I'.ri 1 ) = (i -§’)R<o,n') + £' R(i,n'> + (i -r}')R(£'.o)+ n' R(£' . i) 

- (1 -§')(! - n') R t - (1 -S')n’ R 2 -I 1 v R 3 - V (1 -n')R 4 (3) 

where 


R(£\n') 


x(| ',n') 

y (£'>17') 


is a two -component vector function and £',r 7 * range from 0 to 1 along the con- 
tour of the region. The quantities R^R^.R^R^ are the x, y coordinates of 
the corners of the region and R (0 ,rf), R(l,rj') > R(£’, 0), R(£', 1) are geometric 
functions which describe the shape of the sides of the region (Fig. 2-2). To 
describe the geometry of a region thus requires input of four (x, y) coordinate 
pairs and four sets of functions to describe the sides. Figure 2-4 lists the 
library currently contained in the GIM code. 
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Fig. 2-1 - Flow Domain Example Subdivided into Regions 



Fig. 2-2 - Global and Local Coordinate Systems for a Region 


(C 1) 



1 ( 1 , 0 ) 
( 0 , 0 ) 


Fig. 2-3 - Element Breakdown for a Region 





After determining the coordinates and geometric functions, we can then 
transform to ) space and develop the equations in local coordinates. The 

next step is to subdivide each region into an assemblage of interconnected ele- 
ments as depicted in Fig. 2-3. This is accomplished by simply specifying the 
number of node points in each direction. Each element is developed in its own 
local coordinates(£ , rj) in the same manner as the full domain was transformed 
into regions. All types of stretching transformations are possible with Eq. (3) so 
that irregular grid patterns can be used to enhance resolution. The elements are 
considered isoparametric in that the same functions will be used to approximate 
the flow variables and the element geometry. For example, if linear interpolants 
are used for velocity, pressure, etc., then a typical element shown in Fig. 2-3 
will have linear sides. 

Within each element, any function, Q, will be described by a set of inter- 
polants, S, such that 


Q(l >>i) = s^e.n) Q(i.,n.) (4) 

= S. Q. summed i = 1 , k 
i l 

where Ch is the function at the node points of the element and k is the number 
of nodes on the element. The derivatives of any function Q in (£ ,rj) space is 
then 


dQ i „ dQ i _ , . t , 

W = ~§ F Q i W = W Q i summedl=1 - k < 5 > 

since CT is a point function. 

In order to relate such derivatives back to physical (x, y) space, we use 
the following relation from the calculus. 
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where [ J] is the Jacobian matrix of the transformation, defined as 


[J] 


Ox dy 

w w 

Ox dy_ 
Or) Or] 


(7) 


This is easily obtained by differentiating Eq.(3). The integration over an 
element also requires the area transformation 

dx dy = Aj d^ dr) (8) 

where Aj is the determinant of the Jacobian matrix [ j]. 

We now have all of the geometric information to develop the discretized 
conservation laws for an element. 


2.3 DISCRETIZED EQUATIONS AND ASSEMBLY 

With the coordinates and geometric shape functions now defined for 
each element, we proceed to discretize the differential equations for an in- 
dividual element. The procedure follows that of the Method of Weighted 
Residuals (Ref, 3) in that we require the weighted integral of the differential 
equations to be zero over the domain, R. 


I 


/ w (I? + H + If + H ) = / w£rdA = 0 

R R 


(9) 
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where W is an arbitrary weighting function. Equation (9) must be satisfied 
over the domain R and can be written as a summation of integrals over each 
element 


1=2 J 

6 


(10) 


Equation (10) then represents the point of departure for obtaining the discrete 
conservation equations. The integrals over each element are developed 
by applying the general interpolant shape function, S (Eq. (4)) to the differential 
equation (Eq. (1)). For simplicity we will use the same shape functions for 
each flow variable although this is not a restriction on the approach. 

The dependent variables are thus assumed to be given by the following 
relations over an element 


u = s. u. 

E = S. 

E. 

3 3 

3 

3 

F = S.F. 

H = S. 

H. 

J 3 

3 

3 


( 11 ) 


where summation j = l,k is implied, k 
(k = 4 in Fig. 2-3), and U = |y- • 


is the number of nodes on an element 


For any weight function W^, we put Eq. (11) into Eq. (10) for a single 
element to yield 


■//”■ h°j 


as. 


as. 


+ E. + - 75 -J- F. + S. H. 
9x J 9y J J J 


dx dy 


( 12 ) 


y x 


Now we note that U. E., F., H. are point functions dependent only on time, so 

J J J 3 as. as. 

that we evaluate the integrals in terms of known functions, AAL, S^, , g^ • 

Let the following element matrices be defined. 


13 




ff 


W. S. dxdy 
i 3 y 



■ ff 


W. 

l 



dxdy 
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Equation (12) then becomes 


I 

e 


[A«](U. + H j) + [B i e j ] E + [C*] F 


(summed j = 1, k) 


(13) 


(14) 


Equation (14) is a set of k ordinary differential equations for an ele- 
ment e. The integrals (Eq.(13)) can best be evaluated by invoking the (|,rj) 
transformations (Eqs. (6)— (8)) to yield 

1 1 

[A ij ] = ff W i S j A J d ” 
o o 


1 1 



o o 


The discretized element equations are thus given by Eq.(14) with the 
coefficient matrices given by the geometric integrals Eq.(15). Note that in 
Eq.(14) the j index is summed and the i index counts the node number of 
the equation. 
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Since the choice of W is arbitrary, the integral must be stationary for 
variations in W. The so-called quasi-variational approach determines the 
conditions under which the integral is stationary. That is to say that if 
W' = W + 6W then 


61 



e= 1 


R 


e 


6W Z) dxdy 


0 


Now let 

W = W.. (|, H) U. or 

3 3 


6W = W. 6U. 

3 3 


so that 


61 



W. Z) dxdy 6U. | = 


0 


The for 61 to be stationary, the sum of the coefficients of 6U. must be 
zero. This provides a rationale for the classical assembly procedure of MWR 
derived finite element schemes. 


In mathematics, this can be written as 


where 


I~ a 1 = yi 6 • |a?.i 

L mnj m i L 13 J 


6 . 
n 3 


6 . 
mi 


1 if node i of element e coincides 
with node m of the region R 

0 otherwise 


( 16 ) 


and similarly for the B, C matrices. The rationale for doing this is directly 
derivable from the quasi-variational procedure. 
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A simple example of the assembly procedure is given for the following 
four element domain. 


Node 



Node 1 of the full domain contains only element matrices from element (T) 
while node 2 of the full domain contains a sum of element matrix components 
from elements (T) and (T) . Node 5 of the full domain contains a sum of 
matrix components from all four adjacent elements. 


With the nomenclature of Eq. (15), the full domain equations are 


[A ] 
L mn J 


(U + H ) + [B lE + TC ]F = 0 
' n n L mn 1 n L mn J n 


(17) 


Note again that in Eq. (17), the n index is summed over the number of nodes 
and the m index counts equations for each node. In addition, we must recall 
that U,E,F,H are column vectors themselves. For two-dimensional flow, 
Eq.(17) then represents 4k equations for k nodes in a domain, R. The spatial 
analogs and geometric transformations are all contained in the coefficient 
matrices. Thus' Eq. (17) is a system of ordinary differential equations which, 
with suitable initial data, can be integrated as an initial value problem. 

2 . 4 NOD A LANA LOGS 

At this point, Eq.(17) is general in that the geometric functions are 
arbitrary and the flow variable shape functions and weight functions are 
arbitrary. To this point the development is similar to the general Method 
of Weighted Residuals (MWR) (Ref. 3) with the geometric transformations 
built into the integrals. Readers familiar with MWR will realize that closed 
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form evaluation of the integrals (Eq.(15)) is impossible in general due to 
presence of the Jacobian inverse components for an arbitrary transformation. 

In the GIM code we evaluate all element integrals using Gaussian quadrature 
to maintain generality of geometry. 

If we choose the weight functions W to be identical to the shape functions 
S, then the Eqs. (17) represent the Galerkin approach to finite element nodal 
analogs. Inspection of these equations will reveal that, in general, the spatial 
nodal analogs are implicit in that the unsteady derivatives U are summed over 
all connecting node points. In this work we will refer these ’'finite element?’ 
equations as implicit nodal analogs. Any other classical choice of weight func- 
tions are possible at this point in the analysis, such as collocation, least squares, 
etc . 


We now turn to one of the major developments of the GIM analysis; the 
generation of explicit nodal analogs from the MWR point of departure. In 
this work we will refer to explicit analogs as "finite difference" equations 
not to preclude the fact that implicit finite difference analogs are also possible. 
To consider explicit expressions we turn to Eq. (17) and note that the [A^] 
matrix is the key to rendering the algorithms explicit. If the [A?j] element 
matrices were made diagonal then the [A^.^] full domain matrix will be dia- 
gonal. This would reduce Eq.(l7) to an explicit expression for the unsteady 
derivatives, U. 


To accomplish this goal, we choose the weight functions W to be orthogonal 
to the shape functions, S. 



o o 


Wj S t A j d£ drj = 0 


i/j 


(18) 
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Equations (17) then become 
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where are the diagonal (non- zero) components of 


[A ]. 
l mn J 


This choice of weight functions does not reproduce the ''classical' 1 finite 

element form but rather an explicit nodal analog which we may view as a finite 

difference type of analog. The [B 1 and [C ] matrices now resemble "de- 

' mn j l mn J 

rivative takers" in that is approximated by 
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The specific type of finite difference analog which is produced is a function 
of the element shape functions, S, and the orthogonal weight functions, W. 
To clarify this situation and to show that such expressions do yield proven 
finite difference formulas, the following simple example is given. 


Consider the following four -element, nine -node rectangular region: 
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Ay 
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Full Domain 
Numbering System 


Individual Element 
Numbering System 
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The simplest shape functions are the linear interpolants 


(1-$)(1-U> 

£ (1 -V) 

S= U (1 

Using these shape functions in the orthogonality condition, Eq. (18), we find 
the following weight functions. 


W = 


1 

Ax Ay 


r 0Cj(2/3 -£)(2/3 -»j) > 
a 2 (1/3- f) (2/3 - 7 » 

) 

\ oc 3 (1/3 - 4) (1/3 - r/) 
v a 4 (2/3-4) (1/3 - ^ 


( 20 ) 


where the oc. are arbitrary constants. 

l 1 

Putting these S and W functions into the integral coefficients, Eq. (15), 
and carrying out the integration analytically, we get the following element 
matrices . 



1 

36 


0 0 0 

0 -a^ 0 0 

0 0 0 
0 0 0 -a^ 


( 21 ) 
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36>Ay 


-a x 0 0 tt j ’ 

0 a 2 " a 2 0 

0 -a 3 0 

“4 0 0 -« 4 


( 21 ) 


Now we can assemble the coefficient matrices for the center node (Node 5) 
and examine the resulting analog. Applying the assembly technique, Eq.(16), 
the resulting analog for node 5 is the following: 


( a l a 2 + a 3 a 4^ [ F 5 ^ 5 ] 


( a l~ a 4 ) E 8 + (a 3 +a 4 -a 1 -a 2 )E 5 + (a 2 ~a 3 )E 2 

— 


(a l ~ a Z ) F 6 + (" a i + a z +a 3 ' a 4^ F 5 + ^4 _a 3^ F 4 

— 


= 0 


( 22 ) 


Equation (22) is the general finite difference type of analog resulting from the 
orthogonality condition. Since the ol are arbitrary constants we can choose 
them to give a choice of actual nodal analog. Let 


a. 


= a. 


= 1 


a- 


= a. 


= -1 
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Equation (22) then becomes 


E r -E ? F, -f a 

U 5 +H 5 + -2 S- + -W ± = 0 < 23 > 

Equation (23) is readily recognized as a space centered scheme which 
is unconditionally unstable for an explicit single step time integration. The 
important point is, however, that it is a finite difference scheme. The con- 
clusion is that finite difference schemes (explicit nodal analogs) arise from 
orthogonal weight and shape function sets while classical finite element 
schemes (implicit nodal analogs) result from non-orthogonal choices. 

By selecting the ol in various manners other finite difference schemes 
may be generated. For example let 


a. = 1 + ju; a, = 1 - jLi; a 9 = a . = -1 


MA* < F fe- 2F 5 + F 4> 


then 

. E 8- E 2 . F 6- F 4 . ,x Ax ( E 8- 2E 5 + E 2> . 

5 + 2Ay + 4 Ax 2 4 Ay ‘ 

which is just Eq. (23) with some psuedo -viscous terms added. 


= 0 (24) 


A more pertinent example is a two-step method in which we choose for 

the first step G'=a=a=0;ct = 1 while for the second step a. = = a A 

2341 ^124 

= 0;a^= 1. We then have for the first (prediction step) 


- Ej- F, - F 

Uc + . 5 + - — x 5 = 0 

5 Ax Ay 


(25a) 


while for the second step using provisional values generated above 




Ay 


= 0 


(25b) 
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. Equations (25) are the MacCormack technique (Ref. 6) which has become 
very popular in recent years. Further it should be pointed out that different 
assemblies may be used for different equations yielding combination schemes 
of forward, backward and centered which have been successfully used in prob- 
lems involving meteorological phenomena and convection analyses. 

Some (perhaps not all) finite difference techniques can therefore be gen- 
erated from the finite element point of departure. These schemes are par- 
ticularly attractive for large scale problems where the implicit nature of the 
non-orthogonal weighting functions create an insurmountable obstacle to the 
successful pursuit of the solution. Another important feature is that a great 
deal is known about stability characteristics of some of the explicit schemes. 

An observation which is felt to be of paramount importance is that it is not 
necessary to make the choice of schemes a priori. Rather a general computer 
program can be written in which the user chooses his assembly techniques upon 
input rather than at inception of the coding effort. Thus, a nonorthogonal choice 
may be made if the number of unknowns is small while an orthogonal choice can 
be made if the number of unknowns is large. This is a tendency then toward a 
very general finite difference/finite element fluid mechanics computer code. 

2.5 BOUNDARY CONDITION TREATMENT 

The boundary condition treatment in the GIM code is based on a quasi- 
variational technique. The discussion here will briefly outline the procedure. 

In this discussion we will assume that any number of constraints are 
acting simultaneously. We generalize the constraint statement to read: 

r i = a ik ^i + ^i = 0; k = l > • • • ' K » i = 1, .... I (26) 

There are then K-I independent motions. Rewritting Eq. (26) yields 
K-I K 

X v X *u +p i = ° < 27 > 

k= 1 i =K-I+1 
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where the i motions are considered dependent while the k subscript motions 
are assumed independent. In matrix form 

W|u k } + [O 1^1 + |fj{ = 0 (28) 

where ] is a square (Ixl) matrix. The dependent motions are found ex- 
plicitly from; 

|*i| = - t^] _1 j[L] |u k { + j P || (29) 

The variational statement is 

K| = -W 1 WjsuJ (30) 

After much algebraic manipulation we get the following expression for 
the independent motions. 

U k = [d' 1 ] ju k - [L] [=£ _1 ] (t^' 1 ] p + U k )| (31) 

where 

[D] = [I] + [L] [L] 

the tilde (~) indicates transpose and the hat (a) indicates new values now cor- 
rected for boundary constraints. Equation (31) is used to compute the inde- 
pendent motions. 

As an example, consider the case of inviscid wall tangency 

q.N=nu+nv+nw = 0 (52) 

x y z ' 

where N is the unit normal to the wall. 
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The constraint coefficients are then 


[L] = [n x iy); [=£ ] = y; p = 0 

From Eqs. (29) and (31) we get the following equations for the following equa- 
tions for the unsteady velocity derivatives. 

U k = U k - N k V N - Vi (33) 

i ii ii ' ' 

j =1 

where the N 1 ? are unit normal components at the wall point i (determined from 
geometry). Inspection will show that this maintains flow tangency if the initial 
t = 0 data are tangent to the wall. 

Consider now the multiple constraint situation (internal corner): 



where the flow is constrained to be tangent to two surfaces hence must flow 
parallel to the corner. The equations of constraint are: 


un, -f vn + wn. = 0 

lx ly lz 

un-. + vn, + wn, = 0 

2x Zy Zz 


(34) 


so that 


M = 

[n. 1 

lx ; 00 = 

n ly n izl 


l n 2xl 

in, n, J 
V Zy Zz 


H = 0 


(35) 
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Again Eqs. (29) and (31) give the appropriate expression for corner flow: 

3 

\jY = V] U j (36) 

i l ii ' ' 

j=l 

where i is the nodal point nunber and k = 1, 2, 3 for each velocity component, 
k. 

The are the unit tangent vector components for the corner (determined from 
the geometry). 

The following are the type of nodal point conditions which the GIM code 
can treat. 


1X2* 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 


Condition 

Known or fixed conditions 

Nodes having the same value 
as another 

Stagnation point 

Internal corner flow 

Wall tangency (slip) 

Not currently used 

Not currently used 

Not currently used 

Wall point determined 
from one-sided differences 

Interior flow points. 


The treatment of each type is in terms of the unsteady derivatives, U, and 
are given by the following: 

Type 0 

This is a boundary node at which all flow variables are specified such 
as an upstream inflow to a duct. Mathematically, this condition is 


25 



th = 0 for all nodes i to be specified 

This requires, of course, that the U vector itself be properly initialized to 
the inflow conditions. 


Inflow 
Type 0 



Type 1 

This is a boundary condition for points which are common to more than 
one element such as: 

8 



Nodes 1, 2, 3 and 4 are actually the same point but four values may be deter- 
mined from flow in each of four elements. Mathematically, this is treated 
as follows: 


N 

6 i = Tj £ U i = 1,N (37) 

j=l 


Type 2 


A type 2 node is a stagnation point for which 

U. = 0 for components of U. which 
1 are flow velocities 
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Other components of U such as density and energy are integrated using the 
n-odal analogs. 


Type 2 



Type 3 

Flow in a corner is forced to flow along the corner by imposing condi- 
tions on the components of U which are velocity variables 


^ 3 

Ijk = T^ ^ T 1 ? irf (38) 

l l i i 

j= 1 


where i is the nodal point number and k= 1, 2, 3 for each velocity component. 
1c 

The T^ are the tangent vector components for the corner (determined from 
the geometry and the hat ( A ) indicates corrected values). 


Type 4 



Inviscid flow tangency or free slip conditions are specified by imposing 
the following conditions on the velocity component of U: 


A 



u k - N k V 
1 1 

3=1 


N 1 ? 


ijJ 


(39) 
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where the hL are the unit normal components at the wall point i (determined 
from geometry and the hat ( /v ) indicates corrected values). 



This approach maintains tangency if the initial t- 0 data are tangent to the 
wall. 

Jyp e 8 

Type 8 boundary nodes are those in which flow is allowed to exit the 
domain at an unspecified rate. Such a condition is the downstream wall in 
a supersonic flow where the wall values are determined from upstream 
(backward) differences. Wall points of this type are computed from one- 
sided differences determined by a proper choice of the ol to yield one-sided 
nodal analogs at the boundaries. 



Type 9 

These are interior points in the flow domain which are to be unrestricted 
and determined entirely from the nodal analog of the differential equations. 
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2.6 VISCOUS FORMULATION AND STABILITY 


The full set of differential equations solved by the GIM code is shown in 
Fig. 2-5. These will be referred to as the Navier-Stokes equations, although 
they include continuity and energy conservation as well as momentum. For 
inviscid flow calculations without shock waves, the cr and j terms are set to 
zero to produce the Euler set. Figure 2-6 lists the viscous terms with coeffi- 
cients, /i, X, k for the first and second coefficients of viscosity and the thermal 
conductivity, respectively. These are specified by the user as "laminar" con- 
stants. 

For reasons of numerical stability and capture of strong shock waves, 
additional components of the diffusion coefficients are added automatically by 
the GIM code. Figure 2-7 is a list of currently used Numerical Diffusion 
Cancellation (NDC) coefficients. These are added to the real diffusion coeffi- 
cients. The purpose of these NDC coefficients is to cancel low order trunca- 
tion error terms which arise in the numerics. A forthcoming paper by Spradley 
will present the basic principles of the NDC technique in more detail. 

The differential equations are solved in strong conservation or divergence 
law form. The solution is started at some time t where the entire flow field 
mesh is specified. The unsteady or relaxation of the equations is then done 
using the user-specified nodal analogs. At this time, the pressure is computed 
from the ideal gas law for a single component gas. The code is being modified 
to include multi-gas capability and other additions. 

Stability of the unsteady solution is maintained through addition of the 
viscous cancellation terms. The time step is limited by the classical inviscid 
CFL criterion or the diffusion time criteria, whichever is smaller. Adding 
too much artificial viscosity smears the solution and also reduces the allowable 
time step for stability. The GIM code adds just enough viscosity to balance the 
low order truncation errors. 
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Fig. 2-5 - Three-Dimensional Cartesian Navier-Stokes Conservation Laws 


3 





Fig. 2-7 - Example of the Stability Coefficients in GIM Code 




2.7 THE GIM CODE 


The GIM code is divided into four modules: (1) mesh generation; (2) nodal 
analog assembly; (3) unsteady integration; and (4) data display. The mesh gen- 
eration module accepts boundary geometry data, curve or line formula flags, 
and number of cuts in each coordinate direction. A set of general curvilinear 
coordinate maps is then used to subdivide each region into finite elements. 

Each region which is input is processed and then blended. The output is a set 
of coordinates for each element along with the element coefficient matrices. 

The nodal analog assembly module takes the mesh data from a stored external 
file and performs, via quasi-variational procedure, the assembly of the element 
equations into the full domain equations. At this point, the dynamic storage 
allocation is set up so that the unsteady integration module can integrate with 
virtually unlimited problem size. 

The unsteady integration module performs the actual computation of the 
flow by employing the boundary conditions selected by the user. The nodal 
analog at this point is arbitrary and any one of a number of schemes can be 
selected depending on the problem being analyzed. The solution is marched 
forward in time for a specified number of steps or until a steady state is 
reached. The data display module reads the solution profiles from external 
storage (drum, tape) and prints, plots and maps the flow parameters. Figure 
2-8 is a block diagram illustrating the modular construction of the GIM code. 
Figure 2-9 is a summary of the computer utilization of the code. 

A study is currently underway to determine the feasibility and practi- 
cality of implementing the GIM code on a CDC STAR computer system. The 
pipelining features of such a machine appear to be well suited to the GIM 
formulation since the numerical analogs are written as vector products and 
matrix-vector products. The current STAR features also appear to be best 
for explicit methods rather than implicit schemes. The feasibility study is 
being supported by the NASA-Langley Research Center. 
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I 


Module 1 


Module Z 


Module 3 


Module 4 



Fig. 2-8 - GIM Computer Code Structure 


• OPERATIONAL. ON CDC 7600, UNIVAC 1100 SERIES 
% DYNAMIC STORAGE ALLOCATION 

• 7600 CORE REQUIREMENTS 

• 100Kg SCM 

• 50K g LCM (Nominal) 

• 9 External Files 

• EXECUTION TIMES (7600) 

_4 

• Inviscid 4x10 sec/node - Iter 

• Viscous 6x10^ sec/node - Iter 

• Example - 5000 nodes, 300 Iter 

• Inviscid 10 min CP 

• Viscous 15 min CP 


Fig. 2-9 - Computer Utilization Summary 
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3. DEMONSTRATION PROBLEMS 


Seven sample cases were selected for presentation in this technical 
brief because these cases illustrate the utility of the methods and also reveal 
some interesting results for relatively complex flows. The equations are the 
same for all of the problems, but the geometry and boundary conditions are 
different. The equations are the Euler equations plus certain viscous terms 
for an ideal gas written in strong conservative or divergence form. This 
form permits the use of "shock capturing" algorithms rather than "shock 
fitting" algorithms. The MacCormack scheme is used in the computations 
for illustrative purposes. 

3.1 TWO-DIMENSIONAL TRANSONIC NOZZLE 

This case consists of the single phase two-dimensional transonic flow 
in a converging -diverging nozzle with a finite area inlet and a 15 deg exit using 
air as the fluid. The upstream boundary condition consists of specified Mach 
number, density and pressure along a grid arc. Flow tangency is prescribed 
along the centerline and upper wall. The flow will thus be subsonic upstream 
of the throat, then expand to supersonic conditions in the divergent section of 
the region. Figure 3-1 shows the problem geometry and mesh arrangement 
used in the GIM computation. Quadrilateral elements are used and the method 
selected is the two-step MacCormack scheme. Figure 3-2 shows the steady 
state flow vectors on the physical mesh. The length of the arrow gives relative 
magnitude of the velocity. Figure 3-3 is a computer generated plot of the steady 
state Mach contours in the nozzle. These values have been compared to meas- 
ured profiles (Ref. 7) for this configuration for the sonic line. 
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Fig. 3-1 - Geometry for Transonic Nozzle Case 
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0.6 0.B 1.0 1.2 1.4 1.6 


Fig. 3-3 - Mach Number Contours for Transonic Nozzle Case 
(Comparison with Ref. 7 for Sonic Line) 



3.2 TWO-DIMENSIONAL SHOCK CAPTURING 


This case consists of the supersonic flow in a two-dimensional con- 
verging nozzle section to illustrate the "shock capturing" mode of the method. 
Figure 3-4 shows the mesh and geometry for the problem. The upstream 
boundary condition is a specified Mach number - 4, with tangency being en- 
forced on the centerline and the upper wall. Steady state Mach number con- 
tours are plotted in Fig. 3-5 and compared to a previous shock capturing 
forward marching solution. The MacCormack method was selected for use 
in GIM since the previous shock capturing technique (Ref. 8) used this algo- 
rithm. The figure shows the centerline Mach number versus axial distance 
plotted as the solid line and triangles. The GIM calculations are shown by 
the triangles and captures the shock well with the usual amounts of smearing 
that occur with this approach. The upper wall profile is plotted as dashed 
lines and circles. Comparison with the previous MacCormack solutions is 
seen to be good. The previous SCT solutions were generated via a supersonic 
forward marching code, while the GIM results are generated with an unsteady 
relaxation code. 

3.3 INTERNAL BALLISTICS SIMULATION 

This problem consists of an internal ballistics computation. The real 
propellant burning processes and grain regression were not simulated, but 
rather a model was used with wall mass injection at a specified rate. The 
sketch at the bottom of Fig. 3-6 illustrates the model. A solid wall was used 
as the left-most boundary, with an injection distribution from the "top" boundary 
and flow tangency along the centerline. Upper wall and centerline Mach number 
contours are plotted at steady state versus axial distance. The flow is seen to 
be basically one -dimens ional in the straight section of the region. The wall 
Mach number then rapidly decreases as the conical section, is approached, 
and then expands again up to approximately 1.35. A sharp decrease is then 
experienced as the flow tries to slow down through a relatively flat section 
of the "nozzle." Rapid expansion through the conical section is then the final 
regime of flow. At this time, there are no data or previous computation of 
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Mach Number 



Fig. 3-6 - Two-Dimensional Internal Ballistics Computation Mach 

Number vs Axial Distance (X = .3048 m) 
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this type to be used for comparison, but the solution appears to be physically 
correct. 

3.4 THREE-DIMENSIONAL SUPERSONIC FLOW 

Next we consider three-dimensional inviscid flow in a square nozzle as 

shown in Fig. 3-7. The in-flow Mach number is 2.94 (y = 1.4) at the "entrance" 

of the nozzle. The first 10 X q consists of an expansion region with the final 

10 X being a constant area section. The problem thus consists of a supersonic 
o 

expansion-recompression with two intersecting shock waves as depicted in 
Fig. 3-7. This problem has been solved previously by personnel in the Advanced 
Technology Laboratory (Ref. 9) and at Lockheed-Huntsville. Both of these solu- 
tions utilize a supersonic forward marching technique with a relatively fine mesh 
(21 x 18 x several hundred marching steps). We have computed this problem 
using a two-step MacCormack algorithm with the unsteady GIM code with an 
extremely coarse mesh (11 x 11 x 41). This required 540 time steps to reach 
steady state from a ’’cold" start, i.e., no information about the flow structure. 

The very coarse grid was used for this case to: (1) check out and demon- 
strate the basic three-dimensional GIM code without an excessive amount of 
computer time; and (2) to test the stability of the scheme for a very coarse grid. 
The GIM solution is shown in Fig. 3-8 and compared with the ATL solution. Shown 
are plots of pressure versus axial distance at the upper and lower walls in the 
symmetry plane of the nozzle (see Fig. 3-7). The calculations show excellent 
agreement in the expansion region but deviate considerably in the shockwave 
region. They do, however, exhibit the proper behavior with the quantitative dis- 
agreement due possibly to the very coarse grid in the GIM solution or a compu- 
tation problem with the ATL solution. 

Figure 3-9 is a contour map of pressure in the x-y plane at the outer 
wall (z = 2.0) of the nozzle. The plot shows the shock wave and reflection 
clearly even though it is smeared considerably by the crude grid. This case 
is shown to demonstrate the correctness and applicability of the GIM approach 
to complex fluid dynamic computations. Refinement of the solution is now a 
matter of grid size and computer time. 
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Upper Wall 



(z, y) = 1.5 - 0.5 sin (| + x < 10 
(z, y) = 2 , x >10 


(All units dimensionless) 


Fig. 3-7 - Square Nozzle Three-Dimensional Case 
Expansion-Recompression 



Upper Wall (z = 0) 

■■ ATL solution 

O GIM solution 



Fig. 3-8 - Pressure vs Axial Distance at Upper and Lower Wall of Square 
Nozzle — Comparison with ATL Solution (Ref. 9) 

(X = .3048 m, P = 47.87 N/mZ) 





Fig. 3-9 - Steady Pressure Contours in x-y Plane Showing Shock 

Reflection. (Coarse grid llxllx 41)GIM Code Solution 

(X = Y = .3048 m) 
o o 
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3.5 THREE-DIMENSIONAL, TRANSONIC FLOW 


This problem involves the computation of three-dimensional flow in a 
spinning plug nozzle (SPN). The spinning plug nozzle consists of a tubular 
section which necks down to a smaller tubular section, as shown in Fig. 3-10. 

A plug is inserted into the neck-down region to form a three-dimensional 
convergent-divergent nozzle. The plug is symmetric about its own centerline 
but may be deflected from the centerline of the tubular portion, resulting in a 
nonsymmetric flow field. The flow includes subsonic, transonic, and supersonic 
regimes. 

The flow field was computed with the GIM code using 11,628 nodes. The 
results of the computation are presented in Figs. 3-11 and 3-12. These results 
were used in a contractual study conducted by Lockheed for the U. S. Army 
Missile Research & Development Command (MIRADCOM). 

Figure 3-11 shows the Mach number contours at steady state at the upper 
and lower walls of the configuration as shown in Fig. 3-10. Figure 3-12 is the 
corresponding pressure contours. At present there are no data to compare 
with for this configuration, but the results appear to be physically real. Com- 
parison will be made when the data become available. 

3.6 TWO-DIMENSIONAL SHEAR FLOWS 

The flow field analyzed involved mixing the exhaust from a two-dimensional 
Scramjet afterbody nozzle with freestream. The problem configuration and the 
flow properties of the two flow streams are shown in Fig. 3-13. Both flow 
streams have the same value for the ratio of specific heats (y = 1.27). This 
problem was solved under contract from NASA-Langley Research Center. 

The grid used in GIM computation is shown in Fig. 3-14. Quadrilateral 
elements were used with a two-step MacCormack-type of analog. Solutions 
were obtained at steady state by relaxation of the unsteady equations. Figures 
3-15 and 3-l6 are contour maps of pressure and Mach number, respectively. 
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Fig. 3-14 - Finite Element Grid for Two-Dimensional Configuration 
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at this steady state. Figures 3-17 through 3-19 are plots of the GIM solution 
compared to the Seagull code solution (Ref. 10). Figure 3-17 shows vertical 
pressure distributions through the shear layer at three axial stations down- 
stream of the initial mixing zone. Good agreement is seen except at the upper 
wall for the "close-in" stations. Figure 3-18 is the corresponding Mach number 
plots. 


The upper wall pressure distribution is plotted in Fig. 3-19 and compared 
to the Seagull code and the Lockheed Method of Characteristics (MOC) code. 
Agreement is excellent except near X = 3 X q where the GIM results are apparently 
low. This is due to the GIM code treatment of sharp corners. Seagull and MOC 
use a Prandtl-Meyer type expansion on all sharp corners while the GIM code 
solves only the differential equations. 

Aside from these explainable discrepancies, the GIM code results com- 
pare favorably with the Seagull inviscid solution. The shock, although smeared 
by the GIM code, is located in approximately the same place in both analyses and 
the effects of the shear layer are apparent in the GIM code results. 

3.7 THREE-DIMENSIONAL SHEAR FLOWS 

A three-dimensional shear layer, resulting from the interaction of a 
nozzle exhaust stream with the free stream, both beneath and beside the nozzle, 
was computed using the GIM code. The configuration, shown in Fig. 3-20 con- 
sists of a rectangular nozzle suspended below a body or wing. 

Referring to Fig. 3-20, the nozzle has a sharp 20 deg upward turn at 
X = 0, a sharp 6 deg turn downward at X = X qJ and a sharp 6 deg turn outward 
along the dashed line located on the sidewall. The rectangular nozzle exit 
plane is inclined 30 deg from vertical. The body has a sharp 20 deg turn 
upward at X = X Q « 

The problem was analyzed in three parts: (1) the nozzle; (2) the external 
flow upstream of the nozzle exit, and (3) the entire flow downstream of the exit 
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* plane. The nozzle flow field, which is three-dimensional, was computed using 
the GIM code. It was computed using 2160 nodes and required 570 iterations 
to relax to the steady solution. The external flow was computed using a Prandtl 
Meyer expansion at the 20 deg sharp turn. The nozzle exit plane conditions and 
the freestream expansion were used as input boundary conditions for the down- 
stream shear layer region. 

The shear layer region was analyzed using 7094 nodes and converged in 
260 iterations. A single grid surface is shown in Fig. 3-21 to illustrate the 
grid density and alignment used in the analysis. The results of the analysis 
are presented in Figs. 3-22 through 3-27. A cross section of the problem is 
inserted in each figure to indicate which areas of the flow field the data repre- 
sent. 


Figure 3-22 presents vertical pressure profiles in the plane of symmetry 
for three downstream stations. The profiles illustrate the downward movement 
and dissipation of the shear layer. Figure 3-23 presents pressure profiles 
along the upper wall at constant values of X. The gradients are more severe 
due to the expansion of the free stream, and the 20 deg turn, from 105.8 P Q to 
7.5 P 0 . The shear layer along the upper wall is shown to move outward and 
dissipate. 

Figures 3-24 through 3-26 present pressure contours at three axial 
stations. An outline of the nozzle has been superimposed for reference. 

Figure 3-24 shows the thin shear layer just downstream of the exit plane and 
the freestream pressure gradient due to the Prandtl-Meyer expansion. Figures 
3-25 and 3-26 illustrate how the external flow rolls upward and outward when 
not constrained by the nozzle walls. The dissipation and outward movement 
of the shear layer are also illustrated. Figure 3-27 presents the pressure 
contours on the upper wall, both inside and outside the nozzle. The nozzle 
sidewall with its 6 deg outward turn is apparent in this figure. The freestream 
Prandtl-Meyer expansion appears as a step function. There are no data to com- 
pare this complex three-dimensional solution, but again from physical insight, 
the results appear at least qualitatively correct. 
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Fig. 3-20 - Three-Dimensional Outboard Module Configuration 





Fig. 3-21 
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Three-Dimensional Shear Layer 
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Fig. 



-23 - Lateral Pressure Distribution at Several Axial Stations 
(Z = .3048 m, P = 47.87 N/m 2 ) 










Fig. 3-25 - P ressure Contours at X = 4.6 X Q 

(X = Y = Z = .3048 m, P = 47.87 N/m 2 ) 
' o o o o 
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Fig. 3-26 - Pressure Contours at X = 6.0 X Q 

(X = Y = Z = .3048 m, P = 47.87 N/m 2 ) 
' o o o o 
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Fig. 3-27 - Pressure Contours on Upper Wall , 

(X = Y = Z = .3048 m, P = 47.8 7 N/m ) 
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4. CONCLUDING REMARKS 


4.1 SUMMARY 

The GIM code described in the preceding sections is in production status 
for two- and three-dimensional flows in arbitrary geometric domains. The 
capability of the production code is currently for' laminar flows of nonreacting 
ideal gases. In addition, however, the GIM code is in a continuous development 
state with modifications and updating of capability. The following are among 
the advantages that GIM has over previous approaches. 

• Geometric flexibility of finite element 

• Integral conservation laws 

• Knowledge base of finite difference 

• Boundary conditions via quasi-variational 

• Speed of finite difference 

• No matrix inversion for explicit calculation 

• Unlimited problem size 

• Only methodology which leads to completely 
general fluid mechanics code. 

There are classes of fluid dynamics problems which can be computed 
with much simpler techniques than GIM. However, the flexibility of the GIM 
algorithms make this code applicable to a wide variety of flows. Among the 
classes of problems for which GIM is especially well suited are the following: 

• Flows with subsonic -f supersonic regions 

• Transient analyses 

• Problems requiring full Navier -Stokes equations 

• Complex Geometries 

• Location of three-dimensional shock surfaces. 
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The following summary gives some of the potential areas of fluid dy- 
namics in which the code could be applied, either in its current state or by 
further development. 


4.2 GAS DYNAMICS 


• The code can currently compute a variety of nozzle and exhaust 
plume flow fields for ideal gases in subsonic, transonic, and 
supersonic regimes. 

• Solid rocket motor ballistics computations can be made for 
gas -only flow. Addition of solid particle equations would pro- 
vide the capability of analyzing two phase transonic flow and 
internal gas dynamics for submerged nozzles. 

• Addition of chemistry packages to the GIM code will allow 
computation of reacting gas dynamics for all types of rocket 
motor and plume flow fields. 

• Use of the GIM code for analysis of gasdynamic laser devices 
could prove valuable for three-dimensional flows. This would 
require addition of gas kinetic calculations and a viscous mix- 
ing model. Routines for performing these calculations can be 
added from other present laser codes. 

• The expandable equation sets which can be utilized by the GIM 
code allows inclusion of real viscous terms such that laminar 
or turbulent boundary layer analyses can be made. 

• Inclusion of routines to compute laminar and/or turbulent 
viscosity coefficients would enable a variety of shear flows 
to be computed. Interaction of scramjet exhaust flow with 
freestream air can be treated, including a coupled treatment 
of the shear layer. Low-altitude rocket plume flows can be 
computed, including possible base recirculation with chemistry 
effects and mixing. 


4.3 AERODYNAMICS 

• The GIM code can currently compute subsonic, transonic and 
supersonic flows over arbitrary aerodynamic configurations. 

The code can be used for transonic airfoil design calculations 
although it cannot compete speed-wise with the existing 11 fast 
solvers. 11 Extensions of the code will enable it to be used for 
complete airplane flowfield studies. Aircraft engine flows and 
interactions with the freestream could be computed for scram- 
jet analyses, for example. Also, the potential exists for extend- 
ing this methodology to a hyperbolic marching technique in 
supersonic flow regions. 
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Aerodynamic simulations of wind tunnel tests could be made. 
Effects of tunnel walls and boundary layers can be included 
in computing flow over arbitrary shaped models of aircraft 
configurations. Addition of an appropriate set of boundary 
conditions and a viscosity model will be required. 

Missile aerodynamics presents another area of possible 
application of the GIM code. Flows over high angle of attack 
missiles require a three-dimensional analysis. Computation 
of flows with vortices may be possible by extension of the 
current code. Missile staging aerodynamics present another 
area of possible extension to the current GIM code. 
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